#ifndef DPI_ALGO_H
#define DPI_ALGO_H

#include <math.h>
#include <complex.h>

#include "raylib.h"
#include "FreeImage.h"
#include "fftw3.h"

#if defined(__cplusplus)
extern "C" {            // Prevents name mangling of functions
#endif

///////////////////////////////////////////////////////////////////////////////

//symmetric
typedef enum
{
	SYM_ATTR_ODD = 0, //奇函数
	SYM_ATTR_EVEN, //偶函数
	SYM_ATTR_NEIGHTER
}filter_sym_attr;

///////////////////////////////////////////////////////////////////////////////

#define swapB(a, b) {\
	BYTE temp = a; \
	a = b; \
	b = temp; \
}

typedef int (*func_compare)(BYTE x, BYTE y);

extern int ascending(BYTE x, BYTE y);

extern int descending(BYTE x, BYTE y);

/*
 * 从（无序）数组中，找到第 k 个最小的数（或最大的数)
 * 找到中间的数字
 * kthSmallest(arr, 0, arr_size - 1, arr_size/2, ascending);
 * 返回在数组中的索引, 真正的值是 arr[ret]
*/
extern int kthSmallest(BYTE* arr, int l, int r, int k, func_compare pf);

#define kth_get_min(arr, size) \
	kthSmallest(arr, 0, size - 1, 1, ascending)

#define kth_get_max(arr, size) \
	kthSmallest(arr, 0, size - 1, size, ascending)

#define kth_get_med(arr, size) \
	kthSmallest(arr, 0, size - 1, (size+1)/2, ascending)

// pix_len; default is 3，每个像素占几个字节

extern void gm_translation(Image *ori, Image *chg, int pix_len, int x_offset, int y_offset);

/* is_hon = 1, horizaontal, 水平镜像
 * otherwise, vertical，垂足镜像 */
extern void gm_mirror(Image *ori, Image *chg, int pix_len, int is_hon);

extern void gm_transpose(Image *ori, Image *chg, int pix_len);

///////////////////////////////////////////////////////////////////////////////

/*
* 第3章算法：灰度转换
*/

/*标定图像，让灰度值分布均匀 */

// unsigned char data[size]
extern float data_get_avg(unsigned char* data, int size);

// unsigned char data[size]
// get sigma^2
extern float data_get_sigma2(unsigned char* data, int size, float avg);

/* data[size] to out[size]
 * 值归一化到 [0, 255]
 * */
extern void scale_int_data(unsigned int *data, unsigned char *out, int size);
extern void scale_int_datai(int* data, unsigned char* out, int size);
extern void scale_float_data(float* data, unsigned char* out, int size);

//log(1 + data[i]
extern void log_float_data(float *data, unsigned char *out, int size);

///////////////////////////////////////////////////////////////////////////////

/*
 * 图像反转，反转后的数据写入 chg
 * */
extern void img_negative(FIBITMAP* ori, FIBITMAP* chg);

/*
 * 伽马校正，data用于存放中间结果
 * */
extern void img_gamma(FIBITMAP* ori, FIBITMAP* chg, void* data, float gamma);

//匹配算法后，总的概率数目不是 w*h，重新计算
extern unsigned int histo_get_total(unsigned int* histo);

//计算整幅图的灰度的均值
extern unsigned char histo_get_image_avg(FIBITMAP* dib, unsigned int* histo);

//计算整幅图灰度的标准差
extern unsigned int histo_get_image_sigma(FIBITMAP* dib, unsigned char avg);

/*
 * 从直方图计算变换函数sk, PDF，概率密度分布函数
 * wh: 图像的宽与高的乘积
 * in: histo 直方图概率密度分布函概率密度分布函数
 * out: sk， 变换函数
 * */
extern void get_sk_value(unsigned int* histo, unsigned char* sk, int wh);

/*直方图匹配：最后一步，逆变换 G(z)*/
extern void histo_get_new_mapping(unsigned char* gzq, unsigned char* sk);

/* 把变换函数的值应用到图片中 */
extern void histo_apply_mapping(FIBITMAP* ori, FIBITMAP* chg, unsigned char* sk_data);

/*直方图均衡
 * 变换函数存储在 sk_data[256] 中
 * */
extern void image_histogram_equa(FIBITMAP* ori, FIBITMAP* chg, unsigned char* sk_data);

/* 局部直方图 */

/*
 * 局部直方图
 * type = 3, 3 * 3 邻域
 * type = 5, 5 * 5 邻域概率密度分布概率密度分布函数
 * */
extern void image_histogram_local(FIBITMAP* ori, FIBITMAP* chg, int type);

/* 统计直方图
* type: 3 * 3, 5 * 5
 * avg, sigma: 全局均值，标准差概率?率密度分布函数
 *
 ** if(k0 * avg >= lavg && sigma_l>=k1*sigma_g && sigma_l <= k2*sigma_g)
 *		新的灰度值 = e * 原图的灰度值??率密度分布函数
 * avg: 全局图像均值
 * lavg: 局部均值
 * sigma_g: 图像的标准差概率密度分?率密度分布函数
 * sigma_l: 局部标准差
*/
extern void image_histogram_stats(FIBITMAP* ori, FIBITMAP* chg, int type,\
	unsigned char avg, unsigned int sigma,\
	float e, float k0, float k1, float k2);

/*
 * 两级图像
 * 灰度级 < threshold, value = 0
 * 灰度级 >= threshold, value = 255
 * */
extern void image_binary(FIBITMAP* ori, FIBITMAP* chg, unsigned char threshold);


///////////////////////////////////////////////////////////////////////////////
//空间滤波算法

/*
 * 计算相关/卷积和
 * index: index at image[]
 * sum: spatial correlation/convolution
 * k: 1 based index of mask[dim][dim]
 * val: intensify of current pixel
 * new_val: updated value
 * extra: input/output para
 * */
typedef void (*spatial_filter_func)(unsigned int index, int *sum,\
	unsigned char k, int dim, BYTE val, BYTE* new_val, void *extra);

extern void imgage_filter_common(FIBITMAP* ori, FIBITMAP* chg, int dim,\
	spatial_filter_func filter, void* extra);

/*
 * 均值滤波
 * dim = 3, 3 * 3 averaging mask
 * dim = 5, 5 * 5 averaging mask
 */
extern void image_filter_avg(FIBITMAP* ori, FIBITMAP* chg, int dim);

//中值滤波
typedef enum{
	SPSTAT_MEDIAN = 0,
	SPSTAT_MAX,
	SPSTAT_MIN
}sp_filter_stat;

extern void image_filter_stat_common(FIBITMAP* ori, FIBITMAP* chg,
	int dim, sp_filter_stat type);

#define image_filter_median(ori, chg, dim) \
	image_filter_stat_common(ori, chg, dim, SPSTAT_MEDIAN)

#define image_filter_max(ori, chg, dim) \
	image_filter_stat_common(ori, chg, dim, SPSTAT_MAX)

#define image_filter_min(ori, chg, dim) \
	image_filter_stat_common(ori, chg, dim, SPSTAT_MIN)


//3.6 锐化空间滤波器

//laplacian filter

/*
 * type = 1, used laplacian_coff1
 * type = 2, used laplacian_coff1
 * c = -1
 * otherwise c = 1
 * output:
 *     chg = laplacian result
 * data[h][w]: 存放中间计算结果，因为计算的结果可能是负值
 * data could be NULL
 * */
extern void
image_filter_lapaction(FIBITMAP* ori, FIBITMAP* chg, int type, float* data);

//chg = ori + c * laplacian result
extern void image_filter_lapaction_final
(FIBITMAP* ori, FIBITMAP* chg, int type, float* data);

//Gaussian filter
extern void image_filter_gaussian(FIBITMAP* ori, FIBITMAP* chg, int dim,\
	float sigma);

//3.6.4 grad, Sobel operators
extern void image_filter_grad_sobel_ext(FIBITMAP* ori, FIBITMAP* chg,\
	int* gx, int* gy);
#define image_filter_grad_sobel(ori, chg) \
	image_filter_grad_sobel_ext(ori, chg, NULL, NULL)

/*
 * 6.7-3 - 6.7-9
 * calculate the gradient(magnitude and direction) in RGB space
 * result saved in float data[h][w]
 * */
extern void color_image_gradient(FIBITMAP* color, float* data);

//calc RGB gradient, and merge to data
extern void color_image_gradient_merge(FIBITMAP* color, float* data);

/*
 * 空间滤波器 mask
 * char mask[dim][dim]
 * */
extern
void image_filter_mask(FIBITMAP* ori, FIBITMAP* chg, int dim, char* mask, float* data);

//3.8 fuzzy
extern void fuzzy_intensity_hence(FIBITMAP* ori, FIBITMAP* chg);

extern void fuzzy_spatial_filter_hence(FIBITMAP* ori, FIBITMAP* chg);

///////////////////////////////////////////////////////////////////////////////

/*
 * image_filter_scale_show: just for show
 * image_filter_scale: use it
 * */
extern void image_filter_scale_show(FIBITMAP* dib, float* data);

extern void image_filter_scale(FIBITMAP* dib);

/*
 *	dst = ori + c * add
 */
extern void image_add(FIBITMAP* ori, FIBITMAP* add, FIBITMAP* dst, int c);

extern void
mem_add(BYTE* ori, BYTE* add, BYTE* dst, int c, int w, int h);

//dst = ori + sub
#define mem_addc(ori, add, dst) \
	mem_add(ori, add, dst, 1, w, h)

//dst = ori - sub
#define mem_subc(ori, sub, dst) \
	mem_add(ori, sub, dst, -1, w, h)

/*
 *	dst = ori * mul
 */
extern void image_multi(FIBITMAP* ori, FIBITMAP* mul, FIBITMAP* dst);

///////////////////////////////////////////////////////////////////////////////

//写着玩的，本来想自己实现 FFT
extern void dft(unsigned char* input, float* real, float* imag, int N);
extern void test_dft(void);

extern void dft2d(unsigned char* input, float* real, float* imag, int M, int N);
extern void test_dft2(void);

///////////////////////////////////////////////////////////////////////////////

//ch4 filtering in the frequency domain

typedef enum{
	FFT_SHOW_NORMAL = 0, //[0, 255]
	FFT_SHOW_LOG //1 + logF
}fft_show_type;

typedef enum{
	FREQ_FILTER_IDEAL = 0,
	FREQ_FILTER_BTW, //Butterworth
	FREQ_FILTER_GAUSSIAN,
	FREQ_FILTER_LAPlACE //laplacian
}freq_filter_type;

enum recnotch_mode{
	RECNOTCH_MODE_BOTH = 0,
	RECNOTCH_MODE_HOR, //horizontal
	RECNOTCH_MODE_VER, //vertical
};

typedef struct{
	int w, h;
	void* plan;
	float* image_data; //plan manage it
	fftwf_complex* fft_out; //plan manage it

	BYTE*  image; //image[h][w]
	float* spectrum; //s[h*w]: spectrum data
	BYTE* s_norm; //s_norm[h*w]: normalized data of spectrum data
}fft_data_struc;

/*
 * usage:
 *	  float *image;
 *    fftwf_complex *out;
 *    void* p = fft_init(w, h, &image, &out);
 *
 *    Initialize input signal: fill image[]
 *    FFT: fft_r2c(p, image, out);
 *    IFFT: ifft_c2r(p, out, image);
 *
 *    fft_exit(p);
 *
 *  注意：
 *	  只要宽度和高度不变，就可以无限的用这个 plan
 *      fft_init返回的指针p，只和宽度和高度有关。
 *      所有 fft 处理的数据，用以下4个函数来申请：
 *      fftwf_alloc_real, fftwf_alloc_complex
 *      fftwf_malloc, fftwf_free
 *      这些函数分配的内存保证是对齐的
 *
 *    如果宽度或高度变了，则可以通过调用
 *    fft_init 产生新的计划，只需保留计划指针即可。
 *
 *  in: 图片的宽，高
 *  out:
 *		float image_data[h][w]
 *		fftwf_complex out_data[h][w/2+1]
 *
 *	return:
 *		指向计划的指针
 * */
/*
 * fftwf_free(): release FFT data pointer
 * */
extern void* fft_init(int w, int h, float** image_data, fftwf_complex** out_data);

extern void fft_exit(void* p);

//2D FFT: real to complex
extern void fft_r2c(void* p, float *in, fftwf_complex *out);

//2D IFFT: complex to real
extern void ifft_c2r(void* p, fftwf_complex *in, float *out);

/*
 * 2D FFT: c2c
 * h, w: image size
 * in[h][w]
 * out[h][w]
 */
extern
void fft_c2c(fftwf_complex* in, fftwf_complex* out, int h, int w);

/*
 * 2D IFFT: c2c
 * h, w: image size
 * in[h][w]
 * out[h][w]
 */
extern
void ifft_c2c(fftwf_complex* in, fftwf_complex* out, int h, int w);

extern fft_data_struc* init_fft_data(int w, int h);
extern void destroy_fft_data(fft_data_struc* ptr);


/*
 * 把 FFT 的输出结果输出到完整的复数数组中。 fftwf_complex res[h][w]
 * in: FFT 结果，fftwf_complex* data[h][w/2+1]
 * out: fftwf_complex res[h][w]
 *
 * e.g.
 *  result = malloc(w * h * sizeof(fftwf_complex));
 *  fill_fft_result(data, result, w, h);
 * */
extern void fill_fft_result(fftwf_complex* data, fftwf_complex* res, int w, int h);

/* 计算频谱
 * in: FFT 结果，fftwf_complex* data[h][w/2+1]
 * 把FFT的输出结果，计算频谱后，结果存在 s[]
 * float s[h][w]
 * */
extern void calc_fft_spectrum(fftwf_complex* data, float *s, int w, int h);

/*
 * in:
 *	dib: image handle
 *	center:
 *		1: move spectrum center to [w/2, h/2]
 *		0: keep original
 *	pos: right window pos from 1 to 4
 *  ts: threshold
 *		if ts = 0, ignore it
 *	nw, nh: show as nw * nh
 *
 *	out:
 *	 use destroy_fft_data() to free it
 *	 * */
extern fft_data_struc*
draw_spectrum(FIBITMAP* dib, int center, int pos, int ts, int nw, int nh);

/*
 * in:
 *	data at p->image_data
 *	center:
 *		1: move spectrum center to [w/2, h/2]
 *		0: keep original
 *	pos: right window pos from 1 to 4
 *  ts: threshold
 *		if ts = 0, ignore it
 *	nw, nh: show as nw * nh
 *
 *	out:
 *	 use destroy_fft_data() to free it
 *	 * */
extern void calc_draw_spectrum(fft_data_struc* p, 
	int center, int pos, int ts, int nw, int nh);

#define draw_scaled_spectrum(p, pos, nw, nh) \
	calc_draw_spectrum(p, 1, pos, 0, nw, nh)

/* filter image with BandReject filter
 * r: cutoff frequency
 * wb: width of the band
 * order: Butterworth order
 * reverse = 0, br
 * reverse = 1, bp
 *
 * out:
 *	pfft->image_data: original data
 *	pfft->image: scaled data
 */
extern void br_filter_image(fft_data_struc* pfft, FIBITMAP* dib,
	freq_filter_type type, float r, float wb, int order, int reverse);

/*
 * filter image of dib, result at pfft->image
 * this funciton call fftwf_free(filter) to release filter
 * */
extern void
common_filter_image(fft_data_struc* pfft, FIBITMAP* dib, fftwf_complex* filter);


/* 计算相位：phase angle
 * in: FFT 结果，fftwf_complex* data[h][w/2+1]
 * out: float angles[h][w]
 * */
extern
void calc_fft_angle(fftwf_complex* data, float *angles, int w, int h);

/*
 * 从相角数据重建一副图
 * 频谱 = 1
 * in:
 *  p: reuse the fft plan
 *	float angles[h][w]，图片的相角数据, calc_fft_angle 的计算结果
 *	w, h: image width, height
 *
 * out: float res[h][w] 图片
 * */
extern
void fft_build_from_angle(void *p, float *angles, float *res, int w, int h);


/*
 *  * P157, fomula 4.6-26
 * spatial domain: p >= A + B - 1
 * 为了方便FFT计算，p 还需要是偶数 even
 * Avoid wraparound error
 * 可以在空间域避免缠绕错误
 *
 * 在频率域，这样做可以消除振铃错误
 * Avoid ringing error
 *
 * in: 滤波器的宽高mask_w，mask_h, 图片的宽高 w, h
 * out: 新的宽高
 * */
extern void filter_freq_newsize(int mask_w, int mask_h, int w, int h,\
	int* P, int* Q);

/*
 * 空间滤波器 spatial 变换到频域滤波器 out
 * in:
 *	原始的空间滤波器 spatial[h][w]
 *  频域滤波器的宽P, 高Q
 *	原始空间滤波器的奇偶性 attr
 *
 * 返回值: 频率域滤波器
 * 用完之后，由用户调用 fftwf_free 释放内存
 * */
extern
fftwf_complex* spatial_filter_to_freq(char* spatial, int w, int h,
	int P, int Q, filter_sym_attr attr);

/*
 * 在频率域中滤波
 * image[h][w]: 要滤波的图像
 * filter[Q][P]: 频域滤波器, P, Q：滤波器的尺寸
 *
 * out[h][w]：滤波结果放在这里
 * */
extern void filter_freq_domain(unsigned char* image, int w, int h,
	fftwf_complex* filter, int P, int Q, float* out);

#define filter_freq_nopadding(image, w, h, filter, out) \
	filter_freq_domain(image, w, h, filter, w, h, out)

//b = a * b
extern
void cmulti(fftwf_complex* a, fftwf_complex* b, int w, int h);

/*
 * Ideal LowPass Filter
 * r: D0, cutoff frequency
 * p = w, q = h
 * n: Butterworth order. For other filter, ignore it, default is 2
 *
 * out: filter at frequency domain
 * user call fftwf_free(out) to free resources after finished filtering
 */
extern
fftwf_complex* lpfilter(float r, int p, int q, freq_filter_type type, int n);

/*
 * HighPass Filter
 * r: D0, cutoff frequency
 * p = w, q = h
 * n: Butterworth order
 *
 * out: filter at frequency domain
 * user call fftwf_free(out) to free resources after finished filtering
 */
extern
fftwf_complex* hpfilter(float r, int p, int q, freq_filter_type type, int n);

//filter[] = k1 + k2 * filter[] 
extern void 
hpfilter_highboost(fftwf_complex* filter, float k1, float k2, int p, int q);

/*
 * 4.10 BandReject filter, BR
 * r: D0, cutoff frequency
 * w: Width of band
 * p, q: image width and height
 * n: Butterworth order
 *
 * out: filter at frequency domain
 * user call fftwf_free(out) to free resources after finished filtering
 */
extern fftwf_complex*
brfilter(float r, float w, int p, int q, freq_filter_type type, int n);

extern fftwf_complex*
bpfilter(float r, float w, int p, int q, freq_filter_type type, int n);

/* 4.10 Notch filter
 * Notice: no fill with image and filter
 * 交互过程描述如下：
 *  先计算得到原图像的频谱
 *  如果观察到频谱有关于中轴线对称的K个中心亮点
 *  测量出来亮点位置，就是 notch_center
 *  用这个过滤器乘以原图像的FFT，得到的结果显示出来
 *  观察是否在频谱对应的位置上有相应的黑点（过滤了）
 *
 *	这个过程中，原图像和滤波器都不需要填充，就用原图像的大小
 *	最终得到频域滤波器，用这个滤波器滤波时，需要填充
 *
 * */
typedef struct {
	float x;
	float y;
}ncenter;

typedef struct{
	int k; //center pairs nmuber
	ncenter *c; //c[k]
}notch_center;

/*
 * 之前的坐标，是根据原图测量的
 * 对图像滤波时，原图放大2倍，所以
 * 相应坐标需要放大2倍
 *
 * in:
 *   center: 原图上测量的两点中心坐标，相对于原图左上角
 *   scale: default is 2.0
 *
 * out:
 *   新坐标
 *   用完后需要调用 free_center() 释放内存
 * */
extern notch_center*
change_center(notch_center* center, float scale);

extern void free_center(notch_center* center);

/*
 * 4.10 Notch Reject filter, NR
 * r: D0, cutoff frequency
 * c: center pairs
 * p, q: image width and height
 * n: Butterworth order
 * full_size: if full size filter
 *
 * out: filter at frequency domain
 * user call fftwf_free(out) to free resources after finished filtering
 */
extern
fftwf_complex* nrfilter(float r, notch_center* c, int p, int q,
	freq_filter_type type, int n, int full_size);


extern fftwf_complex*
npfilter(float r, notch_center* c, int p, int q,
	freq_filter_type type, int n, int full_size);

/*
 * Generates rectangular notch (axes) filters.
%	H = recnotch(type, mode, P, Q, W, SV, SH) generates an P-by-Q
%	notch filter consisting of symmetric pairs of rectangles of
%	width W placed on the vertical and horizontal axes of the
%	(centered) frequency rectangle. The vertical rectangles start at
%	+SV and -SV on the vertical axis and extend to both ends of the
%	axis. Horizontal rectangles similarly start at +SH and -SH and
%	extend to both ends of the axis. These values are with respect
%	to the origin of the axes of the centered frequency rectangle.
%	For example, specifying SV = 50 creates a rectangle of width W
%	that starts 50 pixels above the center of the vertical axis and
%	extends up to the first row of the filter. A similar rectangle
%	is created starting 50 pixels below the center and extending to
%	the last row. W must be an odd number to preserve the symmetry
%	of the filtered Fourier transform.
%
%       Valid values of type are:
%
%           1: 'reject'    Notchreject filter.
%
%           0: 'pass'      Notchpass filter.
 * */
extern fftwf_complex* recnorch(int type, enum recnotch_mode mode, int p, int q,
	int w, int sv, int sh);

#define recnorch_default(type, mode, p, q) \
	recnorch(type, mode, p, q, 0, 0, 0)

//filer[h][w] = filter * data
extern
void fft_rrmulti(fftwf_complex* filter, unsigned char* data, float* s,
	int w, int h);

///////////////////////////////////////////////////////////////////////////////
//chapter 5

//Geometric mean filter
extern
void image_filter_geo_mean(FIBITMAP* ori, FIBITMAP* chg, int dim);

/*
 * 5.3.3 local noise reduction filter
 * */
extern void
spatial_local_adaptive(FIBITMAP* ori, FIBITMAP* chg, int dim, float variance);

/*
 * 5.3.3 Adaptive, median filter
 * */
extern void
spatial_median_adaptive(FIBITMAP* ori, FIBITMAP* chg, int smax);

//5.6.3 Estimation by Modeling, example 5.10
extern
fftwf_complex* freq_modeling_a(int w, int h, float a, float b, float T);

//example 5.13, 1/H
extern
fftwf_complex* freq_modeling_b(int w, int h, float a, float b, float T);

//H to 1/H
extern void reciprocal_convert(fftwf_complex *in, int w, int h);

//example 5.11
//calc 1/H
extern
fftwf_complex* freq_degradation_a(int w, int h, float k, float r);

//example 5.12, Winer filter
//rH is reciprocal of H (1/H)
extern fftwf_complex*
freq_degradation_wiener(int w, int h, fftwf_complex* rH, float K);


///////////////////////////////////////////////////////////////////////////////
//ch9

typedef enum
{
	//sqaure
	SE_TYPE_SQUARE = 0, //symmetric, r * r, all is 255

	//line
	SE_TYPE_HOR, //horizontal, dim * 1
	SE_TYPE_VER, //vertical, 1 * dim

	//flat, gray type
	SE_TYPE_DISK //disk
}SE_TYPE;

//binary image, 4-connectivity, 8-connectivity

//4-connectivity, r = 1
#define SE_TYPE_FOUR_CON SE_TYPE_DISK

//8-connectivity, r = 3
#define SE_TYPE_EIGHT_CON SE_TYPE_SQUARE

/* A - B
 * Eroding the image with a square structuring element
 * of size dim * dim filter
 * (default components are all 1 of filter)
 * and saved result to chg
 *
 * bc: background color
 * */
extern void
erosion_mem_common(BYTE* ori, BYTE* chg, int w, int h, SE_TYPE type, int r);

#define erosion_ver(ori, chg, r) \
	erosion_mem_common(ori, chg, w, h, SE_TYPE_VER, r)


#define erosion_mem(ori, chg, w, h, r) \
	erosion_mem_common(ori, chg, w, h, SE_TYPE_SQUARE, r)

//simplify
#define erosion_mems(ori, chg, r) \
	erosion_mem_common(ori, chg, w, h, SE_TYPE_SQUARE, r)

#define erosion_eightcon(ori, chg) \
	erosion_mem_common(ori, chg, w, h, SE_TYPE_EIGHT_CON, 3)

extern void
erosion_image_common(FIBITMAP* ori, FIBITMAP* chg, SE_TYPE type, int r);

#define erosion_image(ori, chg, r) \
	erosion_image_common(ori, chg, SE_TYPE_SQUARE, r)

/* A + B
 * Dilation the image with a square structuring element
 * of size dim * dim filter
 * Image and filer: onverlap by at least one element
 * and saved result to chg
 * */
extern void 
dilation_mem_common(BYTE* ori, BYTE* chg, int w, int h, SE_TYPE type, int r);

#define dilation_mem(ori, chg) \
	dilation_mem_common(ori, chg, w, h, SE_TYPE_SQUARE, 3)

#define dilation_eightcon(ori, chg) \
	dilation_mem_common(ori, chg, w, h, SE_TYPE_EIGHT_CON, 3)

#define dilation_fourcon(ori, chg) \
	dilation_mem_common(ori, chg, w, h, SE_TYPE_FOUR_CON, 1)

extern void
dilation_image_common(FIBITMAP* ori, FIBITMAP* chg, SE_TYPE type, int r);

#define dilation_image(ori, chg, type, r) \
	dilation_image_common(ori, chg, type, r)

extern void morph_open_common(BYTE* ori, BYTE* chg, int w, int h,
	SE_TYPE type, int r);

#define morph_open(ori, chg, type, r) \
	morph_open_common(ori, chg, w, h, type, r)

extern void morph_close_common(BYTE* ori, BYTE* chg, int w, int h,
	SE_TYPE type, int r);

#define morph_close(ori, chg, type, r) \
	morph_close_common(ori, chg, w, h, type, r)

extern void
morph_boundary_extraction(FIBITMAP* ori, FIBITMAP* chg, SE_TYPE type, int r);

#define MORPH_START_LABEL 1
#define MORPH_LAST_LABEL 254
#define MORPH_END_LABEL 255

typedef struct{
	BYTE label; //[1, 254], so max 254 lables
	int pixels; //pixel number of this label(connected component)

	int k0x, k0y; //init point of hole
	//rectangle A, start/end x, y
	int sx, sy;
	int ex, ey;
}morph_label;

/*9.5.2
 * fill a hole of the original image
 * hole position: specified by 'hole' parameter
 *
 * input:
 *  h, w: original image height and width
 *	ori[h][w]: data of original image
 * hole: hole information
 *
 * output:
 *  chg[h][w]
 */
extern void morph_hole_filling(BYTE* ori, BYTE* chg, int w, int h,
	morph_label* hole, SE_TYPE type, int r, BYTE* label_ary);

#define morph_label_background(ori, chg, hole, ary) \
	morph_hole_filling(ori, chg, w, h, hole, SE_TYPE_EIGHT_CON, 3, ary)

#define morph_hole_filling_four(ori, chg, hole) \
	morph_hole_filling(ori, chg, w, h, hole, SE_TYPE_FOUR_CON, 1, NULL)


//9.5.3 Extraction of Connected Components

/*
 * label one connected components
 * morph_label* label: it is in-out parameter
 *
 * input:
 *  h, w: image height and weight
 *  ori[h][w]: image data
 *
 *  label->k0x, label->k0y: coordinate of start point, and
 *    color = ori[k0y][k0x], it is foreground color
 *  label->label: label of this connected components
 *
 * output:
 *  label_ary[h][w]: labeled map
 *
 *  label->pixels: pixel number in this connected components
 *  label->k0x, k0y: init point of hole(with background color)
 *	label->sx, sy, ex, ey: rectangle range of this connected components
 *
 * return value:
 *	k steps
 * */
extern int morph_label_fore_con(BYTE* ori, BYTE* label_ary, int w, int h,
	morph_label* label);


/*
 * label all connected components of foreground
 *
 * input:
 *  h, w: image height and weight
 *  ori[h][w]: image data
 *
 *output:
 *  label_ary[h][w]: label map
 *  morph_label label[MORPH_END_LABEL]: fill it
 *
 *return value:
 *	number of connected components
 * */
extern int morph_label_foreground(BYTE* ori, BYTE* label_ary, int w, int h,
	morph_label* label);

/* Morphological Reconstruction
 * w, h: image width and height
 * dilation = 1: dilation
 * otherwise: erosion
 *
 * output: resultR
 */
extern void morph_reconstruction(BYTE* markerF, BYTE* maskG, BYTE* resultR,
	int w, int h, int dilation);

#define morph_recon_dilation(marker, mask, R) \
	morph_reconstruction(marker, mask, R, w, h, 1)

/*
 * set border value according to boder_invert
 * fill hole: boder_invert = 1, border = 255 - border
 * border clear: boder_invert = 0, border value not change
 * not border: value = 0
 * */
extern void
morph_get_marker(BYTE* ori, BYTE* chg, int w, int h, int boder_invert);

#define morph_get_hole_marker(ori, chg) \
	morph_get_marker(ori, chg, w, h, 1)

#define morph_get_clear_marker(ori, chg) \
	morph_get_marker(ori, chg, w, h, 0)


//invert
extern void complement_area(BYTE* buf, int size, int fill);
#define complement_mem(buf, size) \
	complement_area(buf, size, 255)

///////////////////////////////////////////////////////////////////////////////

//gray-scale morphology

/* A - B
 * w, h: image width and height
 * r: radius of SE
 * ori[h][w] to chg[h][w]
 * */
extern void
erosion_mem_gray_common(BYTE* ori, BYTE* chg, int w, int h, SE_TYPE type, int r);

#define erosion_mem_gray(ori, chg, r) \
	erosion_mem_gray_common(ori, chg, w, h, SE_TYPE_DISK, r)

#define erosion_gray_eightcon(ori, chg) \
	erosion_mem_gray_common(ori, chg, w, h, SE_TYPE_EIGHT_CON, 3)

#define erosion_gray_hor(ori, chg, r) \
	erosion_mem_gray_common(ori, chg, w, h, SE_TYPE_HOR, r)


//A + B
extern void
dilation_mem_gray_common(BYTE* ori, BYTE* chg, int w, int h, SE_TYPE type, int r);

#define dilation_mem_gray(ori, chg, r) \
	dilation_mem_gray_common(ori, chg, w, h, SE_TYPE_DISK, r)

#define dilation_gray_eightcon(ori, chg) \
	dilation_mem_gray_common(ori, chg, w, h, SE_TYPE_EIGHT_CON, 3)

#define dilation_gray_hor(ori, chg, r) \
	dilation_mem_gray_common(ori, chg, w, h, SE_TYPE_HOR, r)

extern void
morph_open_common_gray(BYTE* ori, BYTE* chg, int w, int h, SE_TYPE type, int r);

#define morph_open_gray(ori, chg, r) \
	morph_open_common_gray(ori, chg, w, h, SE_TYPE_DISK, r)

extern void
morph_close_common_gray(BYTE* ori, BYTE* chg, int w, int h, SE_TYPE type, int r);

#define morph_close_gray(ori, chg, r) \
	morph_close_common_gray(ori, chg, w, h, SE_TYPE_DISK, r)

extern void
morph_smooth_gray(BYTE* ori, BYTE* chg, int w, int h, int r);

extern void
morph_gradient(BYTE* ori, BYTE* chg, int w, int h, int r);

//top-hat transformation
extern void
morph_tophat(BYTE* ori, BYTE* chg, int w, int h, int r);

//bottom-hat transformation
extern void
morph_bothat(BYTE* ori, BYTE* chg, int w, int h, int r);

//Morphological Reconstruction
// dilation = 1: dilation
// otherwise: erosion
extern void morph_recon_gray(BYTE* markerF, BYTE* maskG, BYTE* resultR,
	int w, int h, int dilation);

//min = min(a,b)
extern void image_get_min(BYTE* a, BYTE* b, BYTE* min, int size);

//max = max(a,b)
extern void image_get_max(BYTE* a, BYTE* b, BYTE* max, int size);

#define morph_recon_dgray(marker, mask, R) \
	morph_recon_gray(marker, mask, R, w, h, 1)


///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////////////


extern void test_2d_r2c(void);
extern void test_2d_c2c(void);

///////////////////////////////////////////////////////////////////////////////



#if defined(__cplusplus)
}
#endif

#endif
